# allows stargazer to produce output for clmm models 
# (see https://github.com/andrewheiss/civil-society-authoritarian-survival)

load("sccbs_recodes.RData")

fake.clm <- function(x) {
  model.fake <- x
  model.fake$nobs <- x$dims$nobs
  model.fake$call[1] <- call("clm")
  model.fake
}

# estimate mobility and expression models for table
mobility.model <- glmer(exit ~ fb_d*fb_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + 
                          pctunemp + pcthisp + pctasian + medhsvl + (1|community), family = binomial(link = "logit"), nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))

mobility.model.hisp <- glmer(exit ~ fb_d*fb_lag + hisp_d*hisp_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + 
                               pctunemp + pcthisp + pctasian + medhsvl + (1|community), 
                             subset = citizen == 1, 
                             family = binomial(link = "logit"))

mobility.model.housing <- glmer(exit ~ fb_d*fb_lag + fb_d*house_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + 
                               pctunemp + pcthisp + pctasian + medhsvl + (1|community), 
                             subset = citizen == 1,  
                             family = binomial(link = "logit"))

expression.model <- clmm(factor(immig) ~ fb_d*fb_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + 
                           pctunemp + pcthisp + pctasian + medhsvl + (1|community),   subset = citizen == 1, link = "logit")

expression.model.hisp <- clmm(factor(immig) ~ fb_d*fb_lag + hisp_d*hisp_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + pctunemp + pcthisp + pctasian + medhsvl + (1|community),  subset = citizen == 1, link = "logit")

expression.model.housing <- clmm(factor(immig) ~ fb_d*fb_lag + fb_d*house_lag + age + education + tenure + ideology + income + gender + own + kids + hispanic + asian + 
 pctunemp + pcthisp + pctasian + medhsvl + (1|community),  
                                subset = citizen == 1)

stargazer(mobility.model, expression.model)

stargazer(mobility.model, fake.clm(expression.model), out = "table1.html", star.cutoffs = c(.20, .10), digits = 2, notes.append = FALSE, single.row = TRUE)

stargazer(mobility.model.hisp, fake.clm(expression.model.hisp), out = "table2.html", star.cutoffs = c(.20, .10), digits = 2, notes.append = FALSE, single.row = TRUE)

stargazer(mobility.model.housing, fake.clm(expression.model.housing), out = "table3.html", star.cutoffs = c(.20, .10), digits = 2, notes.append = FALSE, single.row = TRUE)

# estimate models for graphs

mobility.model <- glmer(exit ~ fb_d.n*fb_lag.n + age + education + tenure + ideology + income + gender + own + pctunemp + pcthisp + pctasian + medhsvl + kids + hispanic + asian + (1|community), subset = citizen == 1,  family = binomial(link = "logit"), nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))

# simulate QOIs and create plot for mobility model results

beta.hat <- fixef(mobility.model)
vcov <- as.matrix(summary(mobility.model)$vcov)
coefsim <- rmvnorm(10000, beta.hat, vcov)
k <- length(beta.hat)
X <- cbind(age,education,tenure,ideology,
           income,gender,own,pctunemp,pcthisp,
           pctasian,medhsvl,kids,hispanic,asian)

fb.qs <- quantile(fb_d.n, c(0, 1), na.rm = TRUE)
lag.qs <- quantile(fb_lag.n, c(.05, .95), na.rm = TRUE)

rep.row<-function(x,n){
  matrix(rep(x,each=n),nrow=n)
}

means0 <- cbind(1, seq(fb.qs[1], 
                       fb.qs[2], length.out = 1000), lag.qs[1],  
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[1] * seq(lag.qs[1], lag.qs[2], length.out = 1000))

results0 <- coefsim %*% t(means0)
results0 <- apply(results0, 1:2, function(x) invlogit(x))

means1 <- cbind(1, seq(fb.qs[1], 
                       fb.qs[2], length.out = 1000), lag.qs[2], 
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                fb.qs[2] * seq(lag.qs[1], lag.qs[2], length.out = 1000))
results1 <- coefsim %*% t(means1)
results1 <- apply(results1, 1:2, function(x) invlogit(x))

fd <- results1 - results0
means <- apply(fd, 2, mean)
sdUpper <- apply(fd, 2, function(x) quantile(x, .95))
sdLower <- apply(fd, 2, function(x) quantile(x, .05))
z <- seq(min(fb_d.n, na.rm = T), max(fb_d.n, na.rm = T), length.out = 1000)

hist_values <- c(0, diff(sapply(1:1000, function(x) as.numeric(sum(z[x] > fb_d.n, na.rm = T)))))

df <- data.frame(z = z * 100, means = means * 100, 
                 sdUpper = sdUpper * 100, sdLower = sdLower * 100) 

df %>% ggplot(aes(x = z, y = means)) + geom_line() + geom_ribbon(aes(ymin=sdLower,ymax=sdUpper),alpha=0.5) + geom_hline(yintercept = 0) + theme_bw() + labs(x = expression(Delta*" Immigrant (1990 - 2000)"), y = "First Difference") + geom_point(aes(z[1], means[1]), size = 10, shape = 21, fill = "white") + geom_point(aes(z[1000], means[1000]), size = 10, shape = 21, fill = "white") + geom_text(aes(z[1], means[1], label = round(means[1], 0), family = "Palatino Linotype")) + geom_text(aes(z[1000], means[1000], label = round(means[1000], 0), family = "Palatino Linotype")) + theme(plot.title = element_text(size=13, family="Palatino Linotype", face = "bold"), axis.title = element_text(size = 12, family = "Palatino Linotype")) + geom_density(aes(x = z, y = hist_values/50), stat = "identity", fill = "black")

ggsave("pb_flight.pdf", device=cairo_pdf, dpi = "retina")

expression.model <- clmm(factor(immig) ~ fb_d.n*fb_lag.n + age + education + tenure + ideology + income + gender + own + pctunemp + pcthisp + pctasian + medhsvl + kids + hispanic + asian + 
                           (1|community), subset = citizen == 1, link = "logit")

# simulate QOIs and create plot for expression model results
#

beta.hat <- coef(expression.model)
vcov <- vcov(expression.model)[-22,-22]
coefsim <- rmvnorm(1000, beta.hat, vcov)

k <- length(beta.hat)
X <- cbind(age,education,tenure,ideology,
           income,gender,own,pctunemp,pcthisp,
           pctasian,medhsvl,kids,hispanic,asian)

means0 <- cbind(seq(fb.qs[1], fb.qs[2], length.out = 1000), 
                lag.qs[1], 
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                lag.qs[1] * seq(fb.qs[1], fb.qs[2], 
                                length.out = 1000))

results0 <- 1 - (coefsim[,4] - coefsim[,5:21] %*% t(means0))
results0 <- apply(results0, 1:2, function(x) invlogit(x))

means1 <- cbind(seq(fb.qs[1], fb.qs[2], length.out = 1000), 
                lag.qs[2], 
                rep.row(colMeans(X, na.rm = TRUE), 1000), 
                lag.qs[2] * seq(fb.qs[1], fb.qs[2], 
                                length.out = 1000))

results1 <- 1 - (coefsim[,4] - coefsim[,5:21] %*% t(means1))
results1 <- apply(results1, 1:2, function(x) invlogit(x))

fd <- results1 - results0
means <- apply(fd, 2, mean)
sdUpper <- apply(fd, 2, function(x) quantile(x, .95))
sdLower <- apply(fd, 2, function(x) quantile(x, .05))
z <- seq(min(fb_d.n, na.rm = T), max(fb_d.n, na.rm = T), length.out = 1000)
hist_values <- c(0, diff(sapply(1:1000, function(x) as.numeric(sum(z[x] > fb_d.n, na.rm = T)))))

df <- data.frame(z = z * 100, means = means * 100, 
                 sdUpper = sdUpper * 100, sdLower = sdLower * 100) 

df %>% ggplot(aes(x = z, y = means)) + geom_line() + geom_ribbon(aes(ymin=sdLower,ymax=sdUpper),alpha=0.5) + geom_hline(yintercept = 0) + theme_bw() + labs(x = expression(Delta*" Immigrant (1990 - 2000)"), y = "First Difference") + geom_point(aes(z[1], means[1]), size = 10, shape = 21, fill = "white") + geom_point(aes(z[1000], means[1000]), size = 10, shape = 21, fill = "white") + geom_text(aes(z[1], means[1], label = round(means[1], 0), family = "Palatino Linotype")) + geom_text(aes(z[1000], means[1000], label = round(means[1000], 0), family = "Palatino Linotype")) + theme(plot.title = element_text(size=13, family="Palatino Linotype", face = "bold"), axis.title = element_text(size = 12, family = "Palatino Linotype")) + geom_density(aes(x = z, y = hist_values/50), stat = "identity", fill = "black")

ggsave("pb_fight.pdf", device=cairo_pdf, dpi = "retina")
